GÉO-VISUALISATION AVEC R

Action nationale de formation
MATE SHS

Packages & données
Exo 1 - Bienvenue à Sète
Exo 2 - Cartographier l’Occitanie
Exo 3 - Cartographier l’Europe
Exo 4 - Cartographier le Monde
BONUS
- Utiliser des données OSM
- Généraliser un fond de carte
BIBLIO
— R infos —

Cette formation de 3 heures 30 porte sur la visualisation de données géographiques sous R. Y sont abordées les traitements SIG de base, la cartographie thématique (cartes en figurés proportionnels, cartes choroplèthes, cartes de typologie, etc) et des cartes reposant sur des techniques plus avancées comme les cartes sur grille ou les cartes de discontinuité.

Sont ici principalement utilisés les packages sf (manipulation de données spatiales) et cartography (cartographie thématique). Chaque exemple propose une chaine de traitement depuis le chargement des données, leur mise en forme, jusqu’à la mobilisation de méthodes adaptées permettant de répondre à des questions spatiales.

Les exemples proposés sont basés sur différentes échelles, du local (région Occitanie) en passant par l’Europe et le Monde.

Packages & données

Environnement de travail

Pour exécuter les programmes proposés par la formation, installez / mettez à jour R (version 3.4.4 minimum), R Studio. Lancez ensuite R studio et installez les packages au moyen des commandes suivantes (prévoir bien 10 minutes)

Les packages

Installation des packages utilisés pendant la formation.

install.packages("sf")
install.packages("cartography")
install.packages("countrycode")
install.packages("readxl")
install.packages("rnaturalearth")
install.packages("cartogram")
install.packages("reshape2")
install.packages("eurostat")

Optionnel

install.packages("rmapshaper")
install.packages("osmdata")

Version des packages utilisés dans ce document

## [1] "sf (version 0.6.3), cartography (version 2.1.2), countrycode (version 1.0.0), readxl (version 1.1.0), rnaturalearth (version 0.1.0), cartogram (version 0.1.0), reshape2 (version 1.4.3), eurostat (version 3.2.9), rmapshaper (version 0.4.0), osmdata (version 0.0.7.1)"

Packages dédiés à visualisation de données géographiques

Le package cartography permet de créer et intégrer des cartes thématiques dans sa chaine de traitement en R. Il permet des représentations cartographiques tels que les cartes en symboles proportionnels, des cartes choroplèthes, des typologies, des cartes de flux ou des cartes de discontinuité. Il offre également des fonctions qui permettent d’améliorer la réalisation de la carte, comme des palettes de couleur, des éléments d’habillage (échelle, flèche du nord, titre, légende…), d’y rattacher des labels ou d’accéder à des APIs cartographiques.
Le package cartogram comme son nom l’indiqie permet de construire des cartogrames continus ou non continus.

Packages dédiés à la manipulation de données (spatiales ou non)

Le package sf est un package qui permet d’importer, gérer et transformer des données géographiques (gestion des projections, opérations SIG).
Le package readxl permet d’importer des fichiers Excel sous R.
Le package reshape2 permet de transformer et mettre en forme des données sous R.
Le package rmapshaper permet de généraliser le contours d’un fond de carte

Packages fournisseurs de données

Le package eurostat permet d’importer des données d’Eurostat, fournisseur de référence pour les territoires européens.
Le package rnaturalearth permet d’importer les fonds de carte vectoriels de référence pour le Monde (niveaux État ou infra-nationaux) ainsi que des éléments d’habillage (graticules, traits de côte, etc.)
Le package osmdata permet d’importer des données d’OpenStreetMap.

Téléchargez les données

Téléchargez les données sur la page github. Décompressez les données sur votre espace de travail.

https://github.com/riatelab/anfdataviz/blob/master/data.zip

Exo 1 - Bienvenue à Sète

Commandes de base

Définissez votre répertoire de travail.

setwd("Votre_repertoire_de_travail")

Consulter le contenu. Et regarder ce qu’il y a dans le répertoire “data”

list.files()

## [1] "data"       "data.zip"   "index.html" "index.Rmd"  "outputs"   
## [6] "pics"       "prg"

list.files("data")

## [1] "download"  "Europe"    "France"    "Hérault"   "Occitanie" "world"

Import d’une couche SIG dans R avec le package sf. Pacakge basé sur GEOS (Geometry Engine Open Source) et GDAL (Geospatial Data Abstraction Library).

library("sf")

## Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3

communes <- st_read(dsn = "data/Hérault/communes.shp", stringsAsFactors = F)

## Reading layer `communes' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Hérault/communes.shp' using driver `ESRI Shapefile'
## Simple feature collection with 343 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

Voir la table atributaire

communes

## Simple feature collection with 343 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## First 10 features:
##    INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1      34029                        BELARGA   Commune simple        415
## 2      34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3      34032                        BEZIERS   Sous-prfecture       9567
## 4      34082                    COMBAILLAUX   Commune simple        902
## 5      34108                     FRONTIGNAN Chef-lieu canton       3984
## 6      34166                      MONTBLANC   Commune simple       2716
## 7      34066                    CAZEVIEILLE   Commune simple       1619
## 8      34296                      SAUSSINES   Commune simple        630
## 9      34269        SAINT-JEAN-DE-MINERVOIS   Commune simple       3277
## 10     34012                     ARGELLIERS   Commune simple       5065
##    NOM_DEPT                       geometry
## 1   HERAULT POLYGON ((742075 6272005, 7...
## 2   HERAULT POLYGON ((770896 6289405, 7...
## 3   HERAULT POLYGON ((718759 6244261, 7...
## 4   HERAULT POLYGON ((761606 6284489, 7...
## 5   HERAULT POLYGON ((758738 6257679, 7...
## 6   HERAULT POLYGON ((728067 6248191, 7...
## 7   HERAULT POLYGON ((762395 6294296, 7...
## 8   HERAULT POLYGON ((784582 6294234, 7...
## 9   HERAULT POLYGON ((686897 6252592, 6...
## 10  HERAULT POLYGON ((756922 6287661, 7...

head(communes)

## Simple feature collection with 6 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##   INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1     34029                        BELARGA   Commune simple        415
## 2     34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3     34032                        BEZIERS   Sous-prfecture       9567
## 4     34082                    COMBAILLAUX   Commune simple        902
## 5     34108                     FRONTIGNAN Chef-lieu canton       3984
## 6     34166                      MONTBLANC   Commune simple       2716
##   NOM_DEPT                       geometry
## 1  HERAULT POLYGON ((742075 6272005, 7...
## 2  HERAULT POLYGON ((770896 6289405, 7...
## 3  HERAULT POLYGON ((718759 6244261, 7...
## 4  HERAULT POLYGON ((761606 6284489, 7...
## 5  HERAULT POLYGON ((758738 6257679, 7...
## 6  HERAULT POLYGON ((728067 6248191, 7...

View(communes)

L’instruction plot permet d’afficher la couche. L’instruction st_geometry permet d’accéder à la variable définissant les géométries.

plot(st_geometry(communes))

Afficher la couche avec des parametres graphiques

plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)

Introduction au SIG avec R / Manipuler les informations spatiales

Nous cherchons à identifier les communes qui ont une partie de leur territoire situé à moins de 20 km du centre de Sete.

Nous commençons par extraire la commune de Sète.

macommune <- "SETE"
monpoly <- communes[communes$NOM_COMM == macommune,]
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)

Extraitre le centroide de la commune de Sete.

moncentre <- st_centroid(x = monpoly)

## Warning in st_centroid.sf(x = monpoly): st_centroid assumes attributes are
## constant over geometries of x

plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)

Calculer une zone tampon

mydist <- 20000

buff <- st_buffer(x = st_geometry(moncentre), dist=mydist)
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)

Récupérer la liste des communes

communes$buff <- st_intersects(st_geometry(communes), st_geometry(buff), sparse = FALSE)
head(communes)

## Simple feature collection with 6 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##   INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1     34029                        BELARGA   Commune simple        415
## 2     34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3     34032                        BEZIERS   Sous-prfecture       9567
## 4     34082                    COMBAILLAUX   Commune simple        902
## 5     34108                     FRONTIGNAN Chef-lieu canton       3984
## 6     34166                      MONTBLANC   Commune simple       2716
##   NOM_DEPT                       geometry  buff
## 1  HERAULT POLYGON ((742075 6272005, 7...  TRUE
## 2  HERAULT POLYGON ((770896 6289405, 7... FALSE
## 3  HERAULT POLYGON ((718759 6244261, 7... FALSE
## 4  HERAULT POLYGON ((761606 6284489, 7... FALSE
## 5  HERAULT POLYGON ((758738 6257679, 7...  TRUE
## 6  HERAULT POLYGON ((728067 6248191, 7... FALSE

Récupérer le nombre de communes à moins de 20 km

nb <- dim(communes[communes$buff == T,])
nb[1]

## [1] 40

Extraction et affichage des communes

communes20km <- communes[communes$buff == TRUE,]
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)

Ajout d’un titre

plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)
montitre <- paste0("Il y a ", nb[1], " communes situées à moins de ", mydist/1000, "km de ",tolower(macommune))
title(montitre)

Export (si besoin)

st_write(obj = communes20km, dsn = "outputs/resultat.shp")

A vous de jouer !

Reproduisez cette procédure avec la commune et la distance de votre choix.

Pour accéder au listing des communes d’Hérault

communes$NOM_COMM

Exo 2 - Cartographier l’Occitanie

Objectifs

Thématique : Réaliser une carte des catégories Socio professionnelles.
Espace d’étude : Occitanie, communes
Visualisation associée : carte en figurés proportionnels (+ cartogrammes de Dorling)
packages utilisés : sf, cartography, cartogram, readxl

Créer un nouveau fichier “Exo2.R”.

Chargement et préparation des données

Pour bien commencer, il est fondamental de jeter un oeil au fichier d’entrée pour comprendre son organisation et savoir comment l’importer sous R.

Ouvrez le fichier data/France/INSEE/pop-act2554-csp-cd-6814.xls Nous nous intéressons ici à l’onglet “COM_2014”.

library("sf")
library("readxl")

Importez correctement cette feuille du tableau Excel dans R. Les commandes de base si dessous permettent de visualiser le nom des colonnes, de filtrer le contenu de la table (par région) et de visualiser finalement le tableau de données dans R.

sheet <- "COM_2014"
data <- data.frame(read_excel("data/France/INSEE/pop-act2554-csp-cd-6814.xlsx", skip = 15, sheet = sheet)) 
colnames(data)

##  [1] "RR"                            "DR"                           
##  [3] "CR"                            "STABLE"                       
##  [5] "DR16"                          "LIBELLE"                      
##  [7] "csx_rec1taxtypac_rec1rpop2014" "csx_rec1taxtypac_rec2rpop2014"
##  [9] "csx_rec2taxtypac_rec1rpop2014" "csx_rec2taxtypac_rec2rpop2014"
## [11] "csx_rec3taxtypac_rec1rpop2014" "csx_rec3taxtypac_rec2rpop2014"
## [13] "csx_rec4taxtypac_rec1rpop2014" "csx_rec4taxtypac_rec2rpop2014"
## [15] "csx_rec5taxtypac_rec1rpop2014" "csx_rec5taxtypac_rec2rpop2014"
## [17] "csx_rec6taxtypac_rec1rpop2014" "csx_rec6taxtypac_rec2rpop2014"

data <- data[data$RR == "76",] # Ne prendre que les communes d'Occitanie (code 76)
View(data) # vOIR Si le fichier a correctement été importé.

Pour notre analyse nous avons besoin de regrouper plusieurs colonnes. Créez de nouveaux champs en aggrégeant emploi + chômage par CSP. Renommez ces champs de façon explicite. Cela facilitera leur traitement par la suite. Agrégez également le code départemental au code communal afin de créer un identifiant (id) compatible avec le fond de carte IGN.

data$id <- paste0(data$DR, data$CR) # Creation du code commune
data[,"agr"]<- round(data[7] + data[8],0) # Actifs + chomeurs par CSP
data[,"art"] <- round(data[9] + data[10],0)
data[,"sup"] <- round(data[11] + data[12],0)
data[,"int"] <- round(data[13] + data[14],0)
data[,"emp"] <- round(data[15] + data[16],0)
data[,"ouv"] <- round(data[17] + data[18],0)
data[,"ouv_chom"] <- round(data[18],0)
data$total <- data$agr + data$art + data$sup + data$int + data$emp + data$ouv

Importez et visualisez le fond de carte (extrait IGN Geofla communes).

communes <- st_read(dsn = "data/Occitanie/occitanie.shp", stringsAsFactors = F)

## Reading layer `occitanie' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Occitanie/occitanie.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4516 features and 17 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 428145.7 ymin: 6137116 xmax: 848019.7 ymax: 6439628
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

plot(st_geometry(communes),col="white", border="black", lwd=0.1)

La Jointure est une étape cruciale qui permet de lier un fond de carte avec des données par un même identifiant. Il faut que les identifiants entre ces deux ressources soient exactement les mêmes.

Lier un fond de carte avec des données :
la jointure

Réalisez la jointure entre données et fond de carte avec l’instruction merge. Cela présuppose en amont que les identifiants géographiques entre ces deux ressources sont les mêmes (ça tombe bien, c’est le cas !). Identifiez bien le nom des champs sur lesquel doit porter la jointure.

communes <- merge(x = communes, y = data, by.x = "INSEE_COM", by.y = "id", all.x=T)

Pour finir, nettoyez le tableau. On ne garde que les colonnes utiles pour la suite…

colnames(communes)

##  [1] "INSEE_COM"                     "ID_GEOFLA"                    
##  [3] "CODE_COM"                      "NOM_COM"                      
##  [5] "STATUT"                        "X_CHF_LIEU"                   
##  [7] "Y_CHF_LIEU"                    "X_CENTROID"                   
##  [9] "Y_CENTROID"                    "Z_MOYEN"                      
## [11] "SUPERFICIE"                    "POPULATION"                   
## [13] "CODE_ARR"                      "CODE_DEPT"                    
## [15] "NOM_DEPT"                      "CODE_REG"                     
## [17] "NOM_REG"                       "RR"                           
## [19] "DR"                            "CR"                           
## [21] "STABLE"                        "DR16"                         
## [23] "LIBELLE"                       "csx_rec1taxtypac_rec1rpop2014"
## [25] "csx_rec1taxtypac_rec2rpop2014" "csx_rec2taxtypac_rec1rpop2014"
## [27] "csx_rec2taxtypac_rec2rpop2014" "csx_rec3taxtypac_rec1rpop2014"
## [29] "csx_rec3taxtypac_rec2rpop2014" "csx_rec4taxtypac_rec1rpop2014"
## [31] "csx_rec4taxtypac_rec2rpop2014" "csx_rec5taxtypac_rec1rpop2014"
## [33] "csx_rec5taxtypac_rec2rpop2014" "csx_rec6taxtypac_rec1rpop2014"
## [35] "csx_rec6taxtypac_rec2rpop2014" "agr"                          
## [37] "art"                           "sup"                          
## [39] "int"                           "emp"                          
## [41] "ouv"                           "ouv_chom"                     
## [43] "total"                         "geometry"

fields <- c("INSEE_COM", "LIBELLE", "agr", "art", "sup", "int", "emp", "ouv","ouv_chom", "total",  "geometry")
communes <- communes[,fields]
colnames(communes) <-  c("id", "name", "agr", "art", "sup", "int", "emp", "ouv","ouv_chom", "total",  "geometry")
head(communes)

## Simple feature collection with 6 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 556068.9 ymin: 6180702 xmax: 612006.7 ymax: 6219884
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##      id          name agr art sup int emp ouv ouv_chom total
## 1 09001 Aigues-Juntes   5   5   5   9   5   5        0    34
## 2 09002  Aigues-Vives   0  19  24  58  73  63       15   237
## 3 09003     Aiguillon   0   0   5  19  44  44       10   112
## 4 09004        Albiès   5  11  11  11  11  16       11    65
## 5 09005          Aleu   4   0   0   4   8   4        4    20
## 6 09006        Alliat   0   0   0  10   5   5        0    20
##                         geometry
## 1 MULTIPOLYGON (((572559.1 62...
## 2 MULTIPOLYGON (((607979.2 62...
## 3 MULTIPOLYGON (((612006.7 62...
## 4 MULTIPOLYGON (((594651 6188...
## 5 MULTIPOLYGON (((557843.5 61...
## 6 MULTIPOLYGON (((584509.3 61...

Visualisation : carte en figuré proportionnel

Le package cartography

Chargez le package cartography.

library("cartography")

Ce package permet de créer et intégrer des cartes thématiques dans sa chaine de traitement en R. Il permet des représentations cartographiques tels que les cartes en symboles proportionnels, des cartes choroplèthes, des typologies, des cartes de flux ou des cartes de discontinuité. Il offre également des fonctions qui permettent d’améliorer la réalisation de la carte, comme des palettes de couleur, des éléments d’habillage (échelle, flèche du nord, titre, légende…), d’y rattacher des labels ou d’accéder à des APIs cartographiques.

Pour utiliser aisément ce package, plusieurs ressources d’intérêts

La documentation du package qui documente toutes les fonctions du package, accessible directement dans R Studio. Pour celà, vous pouvez taper simplement :

?cartography

La cheat sheet de cartography, qui résume les principales fonctions du package de façon synthétique.
La vignette associée au package, qui présente des réalisations issues de ce package, elle aussi accessible directement dans R studio.
La vignette Faire des cartes - Introduction au package sf qui présente intelligemment les packages sf et cartography à partir d’un cas d’utilisation centré sur Mexico.
Le blog R Géomatique, maintenu par l’auteur de cartography qui met à disposition ressources et exemples d’intérêt liés au package et à la représentation d’information spatiale sous R.

Éléments théoriques

Savoir représenter correctement de l’information dans des cartes thématiques, c’est d’abord savoir caractériser la nature statistique des données.

Comment caractériser les données

Jacques Bertin a notamment théorisé dans son ouvrage de référence “Sémiogie Graphique” (1967) les méthodes visuelles qui permettent de retranscrire efficacement une information statistique sous forme cartographique.

Il identifie des variables visuelles (ou rétiniennes) dont l’utilisation dépend de l’implantation (point, ligne, surface) et la nature statistique de l’information à représenter (association, sélection, ordre, quantité).

Bertin, 1967, niveau des variables rétiniennes

Pour les caractères quantitatifs absolus la variable visuelle adaptée est la taille afin de préserver graphiquement les rapports de proportionnalité existant entre les différentes valeurs du tableau de données.

Visualisation !

Cartographiez le nombre d’actifs par commune grâce à la fonction propSymbolsLayer

plot(st_geometry(communes), col="lightblue",border="white", lwd=0.2) # Fond de carte de la région Occitanie
propSymbolsLayer(x = communes, var = "total",  # Carte en figuré proportionnel 
                 symbols = "circle", col =  "red",
                 legend.pos = "right", border = "grey",
                 legend.title.txt = "Les actifs",
                 legend.style = "c")

top <- sort(data.frame(communes)[,"total"], decreasing = TRUE) 
labelLayer(x = communes[data.frame(communes)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) # Nommer les 15 premières villes d'Occitanie en nombre d'actifs
title("Région Occitanie")

Visualisation : carte combinant figuré proportionnel et typologie

Il est possible d’associer à ces figurés proportionnels la représentation d’un caractère qualitatif de différenciation.

La variable visuelle adaptée pour ce type de données est la couleur différencielle.

Exécutez le code ci-dessous, qui propose le calcul d’une typologie sur la surreprésentation d’ouvriers ou de cadres. Cette typologie peut être représentée simplement avec la fonction typoLayer.

Pour le choix des couleurs, vous pouvez utiliser color picker pour obtenir les couleurs que vous désirez au format hexadécimal présenté ci-dessous.

# Typologie (plus d'ouvriers que de cadres ?)
communes$typo <- "Indéterminé"
communes$r <- communes$ouv / communes$sup
communes$r[is.na(communes$r)] <- 0
communes[communes$r > 1.1,"typo"] <- "Plus d'ouvriers que de cadres"
communes[communes$r < 0.91,"typo"] <- "Plus de cadres que d'ouvriers"

# Choix des couleurs
colouv <- "#dd4e44" # tonalité rouge
colsup <- "#63b269" # tonalité verte

# Cartographie
typoLayer(communes, var="typo", 
          legend.values.order = c("Plus d'ouvriers que de cadres", 
                                  "Plus de cadres que d'ouvriers",
                                  "Indéterminé" ),
          col = c(colouv, colsup, 'grey'), 
          border = NA, legend.title.txt = "Type dominant")

On peut aussi décider d’associer cette représentation à la carte précédente en figuré proportionnel grâce à la fonction propSymbolsTypoLayer :

# Combiner avec la carte précédente
plot(st_geometry(communes), col="lightblue",border=NA)
propSymbolsTypoLayer(x = communes, var = "total", var2 = "typo",
                     symbols = "circle",
                     inches = 0.4,
                     col = c(colouv, colsup, 'grey'),
                     border = "white", lwd = 0.1,
                     legend.var2.values.order = c("Plus d'ouvriers que de cadres", 
                                             "Plus de cadres que d'ouvriers",
                                             "Indéterminé" ),
                     legend.var.pos = "right", 
                     legend.var.title.txt = "Nombre d'actifs",
                     legend.var2.title.txt = "Type dominant")

labelLayer(x = communes[data.frame(communes)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) 
title("Cadres et ouvriers en région Occitanie")

Visualisation : cartogramme de Dorling

A cette échelle géographique on remarque que les cercles se surimposent les uns aux autres. Cela ne nuit pas à la lisibilité générale de la carte mais on peut aussi décider d’utiliser un package et une fonction qui permet d’éviter cela.

Pour ce faire, lancez la librairie cartogram et la fonction cartogram_dorling.

library("cartogram")

# cartogram de Dorling
actifs <- cartogram_dorling(communes, "total", k = 1, m_weight = 1, itermax = 20)

plot(st_geometry(communes), col="lightblue",border=NA)

typoLayer(actifs, var="typo", 
          legend.values.order = c("Plus d'ouvriers que de cadres", 
                                  "Plus de cadres que d'ouvriers",
                                  "Indéterminé" ),
          col = c(colouv, colsup, 'grey'), 
          border = "white", lwd=0.2, legend.title.txt = "Type dominant",
          add = T) 
labelLayer(x = actifs[data.frame(actifs)[,"total"] %in% top[1:15],], txt = "name", halo=TRUE, cex = 0.6, col= "#000000", bg = "#FFFFFF50", overlap = FALSE) 
title("Ouvriers et cadres en Occitanie, 2014")

A vous de jouer !

On soite réaliser une carte sur le chômage des ouvriers en Occitanie. Réalisez une carte avec les parametres suivants :

Mode de représentation : cartograme de Dorling combinée + typologie.
Typologie : 2 catégories (communes ou le taux de chômage des ouvriers est supérieur/inférieur à 20 %)
Cercles proportionnels : population ouvrière (actifs + chomeurs)
Nommer les communes où plus de 500 ouvriers sont au chômage

Voici le code à exécuter pour calculer cette typologie (même démarche que plus haut)

communes$typo <- "Indéterminé"
communes$r <- communes$ouv_chom / communes$ouv
communes$r[is.na(communes$r)] <- 0
communes[communes$r >= 0.2,"typo"] <- "Supérieur à 20 %"
communes[communes$r < 0.2,"typo"] <- "Inférieur à 20 %"

Exo 3 - Cartographier l’Europe

Objectifs

Thématique : Réaliser des cartes de densité de population et de PIB par habitant pour les régions européennes.
Espace d’étude : Europes, NUTS3
Visualisation associée : cartes choroplèthes et cartes en grilles
packages utilisés : sf, cartography, eurostat, reshape2

Créez un nouveau fichier “Exo3.R”.

Chargement et préparation des données

Chargez les géométries qui seront utilisées pour visualiser vos résultats. Il s’agit des géométries de référence des territoires européens adaptés à la cartographie thématique proposés par l’UMS RIATE (ref )

# Fond de carte principal
regions <- st_read(dsn = "data/Europe/GREAT/GEOM_EU34_NUTS13/geom_eu34_nuts13.shp", stringsAsFactors = F)

## Reading layer `geom_eu34_nuts13' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/GEOM_EU34_NUTS13/geom_eu34_nuts13.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1480 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2641758 ymin: 1352640 xmax: 7313157 ymax: 5411284
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

# Elements d'habillage
borders <- st_read("data/Europe/GREAT/OTHERS/borders_EU34.shp", stringsAsFactors = F)

## Reading layer `borders_EU34' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/borders_EU34.shp' using driver `ESRI Shapefile'
## Simple feature collection with 96 features and 8 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 2772939 ymin: 1743120 xmax: 5805008 ymax: 5304194
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

countries <- st_read("data/Europe/GREAT/OTHERS/background_countries.shp",stringsAsFactors = F)

## Reading layer `background_countries' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/background_countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 73 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2600302 ymin: 1249109 xmax: 7346435 ymax: 5455936
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

coastlines <- st_read("data/Europe/GREAT/OTHERS/coast_wide.shp", stringsAsFactors = F)

## Reading layer `coast_wide' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/coast_wide.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 2611445 ymin: 1249569 xmax: 7345439 ymax: 5454917
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

cyprus <- st_read("data/Europe/GREAT/OTHERS/cyprus_north_part.shp", stringsAsFactors = F)

## Reading layer `cyprus_north_part' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/cyprus_north_part.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 6380204 ymin: 1646065 xmax: 6523706 ymax: 1761964
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

remote <- st_read("data/Europe/GREAT/OTHERS/remote_territories_wide.shp", stringsAsFactors = F)

## Reading layer `remote_territories_wide' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/remote_territories_wide.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2657347 ymin: 3613414 xmax: 2881832 ymax: 3836004
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

seaboxes <- st_read ("data/Europe/GREAT/OTHERS/remote_area_background_wide.shp",stringsAsFactors = F)

## Reading layer `remote_area_background_wide' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/remote_area_background_wide.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 23 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2657347 ymin: 2816552 xmax: 2881832 ymax: 4632866
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

background <- st_read("data/Europe/GREAT/OTHERS/background.shp", stringsAsFactors = F)

## Reading layer `background' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/GREAT/OTHERS/background.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2600302 ymin: 1249109 xmax: 7346435 ymax: 5455936
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

Il est toujours utile de visualiser les géométries utilisées et consulter la table attributaire avant de lancer les traitements…

# Ajuster les marges
par(mar = c(0.5,0.5,1.5,0.5)) 
# Visualisation des couches géographiques
plot(st_geometry(background), col = "#e0ecff", border = NA, add =FALSE)
plot(st_geometry(countries), col = "#e6e6e6", border = NA, add =TRUE)
plot(st_geometry(seaboxes), col = "#f7fcfe", border = NA, add = TRUE)
plot(st_geometry(regions),  col = "#7ccdea", border = "#f7fcfe", lwd = 0.15, add=T) # Les representations cartographiques porteront sur cette couche 
plot (st_geometry(remote), col = "#e6e6e6", border =NA, add=TRUE)
plot (st_geometry(borders), col = "#f7fcfe",lwd = 0.5, add = TRUE) 
plot (st_geometry(cyprus), col = "#f7fcfe", border = NA, add = TRUE)
plot (st_geometry(coastlines), col = "#d2dbef", lwd = 0.15, add = TRUE)
plot(st_geometry(seaboxes), col = NA, border = "#bbbdc0", lwd = 0.15, add = TRUE)

# Table attributaire de la couche cible pour les traitements cartographiques
head(regions)

## Simple feature collection with 6 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 4655611 ymin: 2655076 xmax: 4855527 ymax: 2823072
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
##      id                    name level0 level1 level2 level3 level23
## 1 AT111        Mittelburgenland      0      0      0      1       0
## 2 AT112          Nordburgenland      0      0      0      1       0
## 3 AT113           Südburgenland      0      0      0      1       0
## 4 AT121 Mostviertel-Eisenwurzen      0      0      0      1       0
## 5 AT122    Niederösterreich-Süd      0      0      0      1       0
## 6 AT123            Sankt Pölten      0      0      0      1       0
##                         geometry
## 1 MULTIPOLYGON (((4819333 274...
## 2 MULTIPOLYGON (((4819333 274...
## 3 MULTIPOLYGON (((4808360 271...
## 4 MULTIPOLYGON (((4729815 281...
## 5 MULTIPOLYGON (((4762433 279...
## 6 MULTIPOLYGON (((4768449 279...

Ce bloc de code est dédié est à la recolte des données disponibles sur le site d’Eurostat grâce au package eurostat. La procédure à suivre pour télécharger ses propres données est la suivante :

Consultez le contenu de la base de données Eurostat et identifiez la table qui contient vos indicateurs d’intérêt (cf ci-dessous)
Importez la librairie eurostat.
Téléchargez la table d’intérêt avec la fonction get_eurostat. C’est un tibble.
Consultez les dimensions du tibble avec la fonction str()
Filtrez les dimensions du tibble afin de réduire sa structure à deux dimensions (objets géographiques * indicateur par exemple).
Transformez le tibble en dataframe avec la fonction dcast(). Dans notre cas de figure nous conservons en ligne les objets géographiques et en colonne la dimension temporelle de l’indicateur d’intérêt.
Quand tous les indicateurs ont été téléchargés et préparés, les joindre à vos géométries de référence (regions dans notre cas de figure).

Base de données régionale d'Eurostat></img></p>
<div class=

library("eurostat")
library("reshape2") 

# POPULATION
var <- "demo_r_pjangrp3" # Table Eurostat d'intérêt
data <- get_eurostat(var, time_format = "num") # Telecharger la table ESTAT
data <- subset(data, data$sex == "T") # Filtre des dimensions du tibble
data <- subset(data, data$age == "TOTAL") # Filtre des dimensions du tibble
data <- dcast(data, geo ~ time, value.var = "values") # Transformation en dataframe des dimensions restantes
colnames(data) <- c("geo","POP_2014","POP_2015","POP_2016","POP_2017") # Renommer les colonnes du dataframe de façon explicite

# SURFACE
var <- "reg_area3"
area <- get_eurostat(var, time_format = "num") # Telecharger la table ESTAT
area <- subset(area, area$landuse == "TOTAL") # Filtre des dimensions de la table
area <- dcast(area, geo ~ time, value.var = "values")
colnames(area) <- c("geo","AREA")

# PIB
var <- "nama_10r_3gdp"
gdp <- get_eurostat(var, time_format = "num")
gdp <- subset(gdp, gdp$unit == "MIO_EUR")
gdp <- dcast(gdp, geo ~ time, value.var = "values")
fields <- c("geo", "2014", "2015", "2016") # On ne garde que les valeurs de PIB correspondant aux valeurs de population disoponibles
gdp <- gdp[,fields]
colnames(gdp) <- c("geo","GDP_2014","GDP_2015","GDP_2016")

# Jointure
regions <- merge(x = regions, y = data, by.x = "id", by.y = "geo", all.x=T) # population
regions <- merge(x = regions, y = area, by.x = "id", by.y = "geo", all.x=T) # surface
regions <- merge(x = regions, y = gdp, by.x = "id", by.y = "geo", all.x=T) # surface

Importez le jeu de données créé suite à cette procédure, associé aux géométries européennes de réference.

regions <- st_read("data/Europe/regions_data.shp", stringsAsFactors = F)

## Reading layer `regions_data' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Europe/regions_data.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1480 features and 15 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2641758 ymin: 1352640 xmax: 7313157 ymax: 5411284
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

Les données sont prêtes pour la visualisation !

Pour choisir l’année de référence sur laquelle portera la visualisation, regardez le nombre de valeurs manquantes…

apply(regions,2,function(x) sum(is.na(x)))

##       id     name   level0   level1   level2   level3  level23 POP_2014 
##        0        0        0        0        0        0        0        2 
## POP_2015 POP_2016 POP_2017     AREA GDP_2014 GDP_2015 GDP_2016 geometry 
##        1        1        1        1      111      113     1161        0

… Et calculez le ratio sur lequel va porter la visualisation en conséquence

regions$DENS_2017 <- regions$POP_2017 / regions$AREA

Visualisation : cartes choroplèthes

La variable visuelle adaptée pour représenter des données quantitatives relatives est la couleur ordonnée (ou niveaux de gris). Le traitement de ce type de données implique au prélable de définir le nombre de classes, une méthode de discrétisation et une palette de couleurs sur laquelle portera la représentation.

La discrétisation consiste à découper une série statistique en classes de valeurs. Il n’existe pas de méthode optimum. Chaque méthode donnera une carte différente, plus ou moins conforme à la distribution de départ. L’agrégation de données en classes introduit de fait une erreur ou une distorsion dans la perception de cette distribution. Le choix d’une méthode de discrétisation dépend donc des propriétés de la distribution mais aussi des objectifs cartographiques que l’on s’est fixés.

Dans cartography, plusieurs méthodes de discrétisation sont proposées. Jetez un oeil à la documentation de la fonction getBreaks pour en savoir plus.

On peut maintenant observer avec un graphique simple les effets de plusieurs méthodes de discrétisation sur notre variable de densité de population.

Compte-tenu du nombre d’unités territoriales NUTS3 en Europe, nous choisissons 10 classes pour guider la discrétisation.

library(cartography)

# 4 graphiques par plot
par(mfrow=c(2,2))

# Methode discretisation (quantiles)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "quantile")
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Quantiles")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)

# Methode discretisation (intervalles égaux)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "equal")
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Intervalles égaux")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)

# Methode discretisation (moyenne-écart type)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "msd", k=1, middle = TRUE)
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Moyenne/écart-type")
med <- median(regions$DENS_2017, na.rm = T)
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)

# Methode discretisation (geom)
breaks <- getBreaks(v = regions$DENS_2017, nclass = 10, method = "geom")
hist(regions$DENS_2017, probability = TRUE, breaks = breaks, col = "#F0D9F9", xlim = c(0,4000), main = "Progression géométrique")
abline(v = median(regions$DENS_2017, na.rm = T), col = "blue", lwd = 3)
abline(v = mean(regions$DENS_2017, na.rm = T), col = "red", lwd = 3)

Explorez la documentation de la fonction carto.pal de cartography pour choisir votre palette. Utilisez l’instruction display.carto.all pour visualiser les palettes disponibles. Nous choisissons ici la palette “red.pal”

# Choix de la palette
display.carto.all(10)

cols <- carto.pal(pal1 = "red.pal", n1 = 10)

Tout est prêt ! Cartographions cet indicateur avec la méthode des quantiles, la palette “red” en 10 classes. Associons des couches d’habillage. Et mettons en page la carte (flèche du nord, échelle, sources) grâce à l’instruction layoutLayer de cartography.

# Ajuster les marges et les paramètes de visu
par(mar = c(0.5,0.5,1.5,0.5)) 
par(mfrow=c(1,1))

# Mise en page (back)
plot(st_geometry(background), col = "#e0ecff", border = NA, add =FALSE)
plot(st_geometry(countries), col = "#c6c4c4", border = NA, add =TRUE)
plot(st_geometry(seaboxes), col = "#f7fcfe", border = NA, add = TRUE)

# Cartographie (ratio)
choroLayer(x = regions, 
           var = "DENS_2017",
           method ="quantile",
           nclass = 10,
           col = cols,
           border = NA,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "Densité de population, 2017 (hab/km²)",
           legend.values.rnd = 0)

# Mise en page (top)
plot (st_geometry(remote), col = "#e6e6e6", border =NA, add=TRUE)
plot (st_geometry(borders), col = "#f7fcfe",lwd = 0.5, add = TRUE) 
plot (st_geometry(cyprus), col = "#f7fcfe", border = NA, add = TRUE)
plot (st_geometry(coastlines), col = "#d2dbef", lwd = 0.15, add = TRUE)
plot(st_geometry(seaboxes), col = NA, border = "#bbbdc0", lwd = 0.15, add = TRUE)

# Sources + elements d'habillage
layoutLayer(title = "Densité de population dans les NUTS3 européens, 2017",
            author = "Réalisation : Serial Mappers, 2018", sources = "Source, Eurostat, 2018",
            scale = 200, tabtitle = TRUE,
            frame = FALSE,theme = "red.pal",
            south = TRUE)

A vous de jouer !

Reproduisez cette carte avec les paramètres suivants :

Discrétisation : progression géométrique
Nombre de classes : 8
Palette de couleurs : Bleue

Visualisation : cartes carroyées

La méthode du carroyage consiste à découper l’espace géographique en un maillage formé de carrés réguliers dans une projection donnée. Elle permet de s’affranchir de l’arbitraire et de l’irrégularité d’un découpage administratif, comme c’est ici le cas avec l’Europe au niveau NUTS3. La donnée est répartie sur ce quadrillage régulier au prorate de la surface représentée.

Comment faire ? L’instruction getGridLayer du package cartography permet de construire ces grilles régulières. Expolorez la documentation pour vous sensibiliser aux arguments de la fonction.

La procédure est la suivante :

Définir le pas de grille (en mètres), son type et les variables qui sont transformées
Pour les calculs de densité, convertir la surface en km².
Pour les autres indicateurs, combiner les données de la façon suivante : ratio = stock1/stock2
Choisir sa palette, son nombre de classes et sa méthode de discrétisation
Représenter la grille avec la fonction choroLayer de cartography.

# Transformation des donnees dans une grille reguliere de 50 km
mygrid50k <- getGridLayer(x = regions, cellsize = 50000 * 50000, 
                       type = "regular", var = c("POP_2017","POP_2014","GDP_2014"))


# Transformation des donnees dans une grille reguliere de 100 km
mygrid100k <- getGridLayer(x = regions, cellsize = 100000 * 100000, 
                          type = "regular", var = c("POP_2017","POP_2014","GDP_2014"))

# Conversion en km²
mygrid50k$densitykm <- mygrid50k$POP_2017 * 1000 * 1000 / mygrid50k$gridarea 
mygrid100k$densitykm <- mygrid100k$POP_2017 * 1000 * 1000 / mygrid100k$gridarea 


# Cartographie
par(mfrow=c(1,2)) # 2 cartes
par(mar = c(0.5,0.5,1.5,0.5)) # Ahuster les marges

cols <- carto.pal(pal1 = "blue.pal", pal2 = "red.pal",  n1 = 5, n2 = 5)  
  
choroLayer(x = mygrid50k, var = "densitykm", 
           border = "grey80",col = cols, method ="quantile",
           nclass = 10, legend.pos = "topright", legend.values.rnd = 1,
           legend.title.txt = "Densité de population, 2017")

layoutLayer(title = "Densité de population dans une grille régulière de 50km",
            author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
            scale = 200,
            frame = TRUE,
            col = "red",
            coltitle = "white",
            south = TRUE)

choroLayer(x = mygrid100k, var = "densitykm", 
           border = "grey80",col = cols, method ="quantile",
           nclass = 10, legend.pos = "topright", legend.values.rnd = 1,
           legend.title.txt = "Densité de population, 2017")

layoutLayer(title = "Densité de population dans une grille régulière de 100km",
            author = "Serial Mappers, 2018", sources = "Eurostat, 2018",
            scale = 200,
            frame = TRUE,
            col = "red",
            coltitle = "white",
            south = TRUE)

A vous de jouer !

Créez une nouvelle carte carroyée avec les paramètres suivants :

Indicateur : PIB / habitant 2014
Résolution de grille : 100 km
Type de grille : hexagones (nb : il faut charger la librairie sp au préalable)
Discrétisation : quantiles (8 classes)
Palette de couleurs : Double palette (vert/rouge)

Exécutez ce code qui permet au prélable de préparer la grille et retirer les valeurs manquantes qui nuiraient au calcul.

library(sp)

# Retirer les valeurs manquantes
mygrid100k <- regions
mygrid100k <- mygrid100k[!is.na(mygrid100k$GDP_2014),]
mygrid100k <- mygrid100k[!is.na(mygrid100k$POP_2014),]

# Creation de la grille
mygrid100k <- getGridLayer(x = mygrid100k, cellsize = 100000 * 100000, 
                          type =  "hexagonal", var = c("POP_2014","GDP_2014"))

# Calcul PIB par habitant
mygrid100k$GDP_HAB <- mygrid100k$GDP_2014 * 1000000 / mygrid100k$POP_2014

Exo 4 - Cartographier le Monde

Objectifs :

Thématique : Les inégalités d’IDH et/ou d’espérance de vie
Espace d’étude : Le Monde, pays.
Visualisation associée : carte chorplèthe, discontinuités
packages utilisés : sf, cartography, rnaturalearth, readxl, countrycode

Import du fond de carte

méthode 1 : à la main (après téléchargement sur le site Natural Earth)

library("sf")
countries <- st_read(dsn = "data/world/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F)

Méthode 2 : directement via R

url <- "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip"
download.file(url = url, destfile = "data/download/countries.zip")
unzip("data/download/countries.zip", exdir = "data/download", unzip = "internal")
file.remove("data/download/countries.zip")
countries <- st_read(dsn = "data/download/ne_110m_admin_0_countries.shp", stringsAsFactors = F)

Méthode 3 : directement via R

library("rnaturalearth")
countries <- ne_countries(scale = 50, type = "countries", continent = NULL,
                          country = NULL, geounit = NULL, sovereignty = NULL,
                          returnclass = "sf")

Le fond de carte peut être visualisé avec l’instruction plot

plot(st_geometry(countries))

On ne conserve que les champs utiles. On les renomme.

countries <- countries[,c("adm0_a3", "admin", "geometry")]
colnames(countries) <- c("id","name","geometry")

On Supprime l’Antarctique

countries <- countries[countries$id != "ATA",]

On supprime les lignes vides

countries <- countries[!is.na(countries$id),]

plot(st_geometry(countries))

Import des données atributaires

Import des données sur l’éspérance de vie (source : Banque mondiale)

library("readxl")
file <-"data/world/worldbank/API_SP.DYN.LE00.IN_DS2_en_excel_v2_10081535.xls"
sheet <- "Data"
lifexp <- data.frame(read_excel(file, skip = 3, sheet = sheet))
lifexp <- lifexp[,c("Country.Code","Country.Name","X2016")]
colnames(lifexp) <- c("id","name","lifexp2016")
lifexp$lifexp2016 <- as.numeric(as.character(lifexp$lifexp2016))

Import des données sur l’Indice de deveoppement humain (source : Human Development Report)

file <-"data/world/hdr/Human Development Index (HDI).csv"
hdi <- read.csv2(file = file, sep = ",",skip = 1)
hdi <- hdi[,c("Country","X2016")]

head(hdi,4)

##        Country X2016
## 1  Afghanistan 0.494
## 2      Albania 0.782
## 3      Algeria 0.753
## 4      Andorra 0.856

Le tableau de données sur l’IDH ne possède pas de codes géographiques, mais uniquement le nom des pays. Cela rend la jointure avec le fond de carte compliquée. Une solution : le package countrycode

library("countrycode")
hdi$id <- countrycode(hdi$Country, 'country.name', 'iso3c')

## Warning in countrycode(hdi$Country, "country.name", "iso3c"): Some values were not matched unambiguously:  Eswatini (Kingdom of)

Un peu de mise en forme. Les codes sont bien là. C’est presque magique… Mais attention, dans la pratique, une vérification manuelle s’impose.

hdi <- hdi[,c("id","Country","X2016")]
colnames(hdi) <- c("id","name","hdi2016")
hdi$hdi2016 <- as.numeric(as.character(hdi$hdi2016))  
head(hdi,4)

##    id         name hdi2016
## 1 AFG  Afghanistan   0.494
## 2 ALB      Albania   0.782
## 3 DZA      Algeria   0.753
## 4 AND      Andorra   0.856

Cartographie

Chargement du package

library("cartography")

Changer la projection du fond de carte (voir le site spatialreference

countries <- st_transform(countries, "+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")
plot(st_geometry(countries))

Jointures des deux fichiers de données

countries <- merge(countries, lifexp, by="id", all.x=TRUE)
countries <- merge(countries, hdi, by="id", all.x=TRUE)
countries <- countries[,c("id","name.x","lifexp2016","hdi2016","geometry")]
colnames(countries)[2] <- "name"

On travaille sur l’indice de développement humain

analyse de la distribution

var <- countries$hdi2016
hist(countries$hdi2016, probability = TRUE, nclass = 30)

On choisit une discretisation par la méthode des quantiles sans classe centrale

breaks <- getBreaks(v = var, nclass = 6, method = "quantile")
cols <- carto.pal(pal1="blue.pal", n1 = 3, pal2 = "orange.pal", n2 = 3, middle = FALSE, transparency = FALSE)

Réalisation de la carte

layoutLayer(extent = countries,
            bg="#EEEEEE",
            title = "Indice de développement humain en 2016",
            scale = NULL,
            sources = "Human development Report, 2018",
            author = "les serial mappers"
            )
choroLayer(x = countries, 
           var = "hdi2016",
           breaks = breaks,
           col = cols,
           border = "white",
           lwd=0.2,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "IDH, 2016",
           legend.values.rnd = 2)

Et si on ajoutait des lignes de discontinuités à cette carte?

L’instruction getBorders (du package cartography) permet d’extraire les frontières entre les unités spatiales.

borders <- getBorders(countries)
plot(st_geometry(borders), col="black", lwd=1)

Quid des discontinuités maritimes ? Par exemple entre l’Italie et la Tunisie ? Ou entre le Maroc et l’Espagne ? L’instruction getOuterBorders (package cartography) permet de les générer. Cette instruction peut prendre un peu de temps.

outer <- getOuterBorders(x = countries, res = 100000, width=500000)
plot(st_geometry(borders), col="#CCCCCC", lwd=1)
plot(st_geometry(outer), col="red", lwd=1, add=T)

On assemble les deux couches avec rbind

b <- rbind(borders,outer)

Tout est prêt pour réaliser une carte de discontinuités

# 1 - la carte choroplèthe
layoutLayer(extent = countries,
            bg="#EEEEEE",
            title = "Indice de développement humain en 2016",
            scale = NULL,
            sources = "Human development Report, 2018",
            author = "les serial mappers"
            )
choroLayer(x = countries, 
           var = "hdi2016",
           breaks = breaks,
           col = cols,
           border = "white",
           lwd=0.2,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "IDH, 2016",
           legend.values.rnd = 2)

# 2 - les discontinuités

discLayer(x = b, df = countries,
          var = "hdi2016", col="black", nclass=4,
          method="quantile", threshold = 0.25, sizemin = 0.2,
          sizemax = 8, type = "abs", 
          legend.title.txt = "Discontinuities\nabsolues",
          legend.pos = "bottomleft", add=TRUE)

BONUS

Utiliser des données OSM

library("sf")
library("cartography")

depts <- st_read(dsn = "data/France/IGN/DEPARTEMENT.shp", stringsAsFactors = F)

## Reading layer `DEPARTEMENT' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/France/IGN/DEPARTEMENT.shp' using driver `ESRI Shapefile'
## Simple feature collection with 96 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 99217.1 ymin: 6049646 xmax: 1242417 ymax: 7110480
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

colnames(depts)

##  [1] "ID_GEOFLA"  "CODE_DEPT"  "NOM_DEPT"   "CODE_CHF"   "NOM_CHF"   
##  [6] "X_CHF_LIEU" "Y_CHF_LIEU" "X_CENTROID" "Y_CENTROID" "CODE_REG"  
## [11] "NOM_REG"    "geometry"

Paris <- depts[depts$CODE_DEPT == 75,]
plot(st_geometry(Paris))

Paname <- getTiles(x = Paris, type = "osm", crop = TRUE)
Paname <- getTiles(x = Paris, type ="opencycle", crop = TRUE)
tilesLayer(Paname)
plot(st_geometry(Paris), lwd=4, add=T)

library("osmdata")

## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright

prj <- st_crs(Paris)

bbox <- st_bbox(st_transform(Paris,4326))

q <- opq(bbox = bbox , timeout = 2000) %>% add_osm_feature(key = 'man_made', value = 'surveillance')
cameras <- osmdata_sf(q)$osm_points
cameras <- st_transform(cameras, prj)

cameras$ok <- st_intersects(st_geometry(cameras), st_geometry(Paris), sparse = FALSE)
cameras <- cameras[cameras$ok == T,]

Paname <- getTiles(x = Paris, type ="cartodark", crop = TRUE)
tilesLayer(Paname)
plot(st_geometry(Paris), lwd=1, border="white", lty=2, add=T)
plot(st_geometry(cameras), pch=20, col="red", add=T)

Généraliser un fond de carte

Bla bla sur la généralisation cartographique

library("rmapshaper")

Import et affichage d’un fond de carte

countries <- st_read(dsn = "data/world/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F)

## Reading layer `ne_110m_admin_0_countries' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/world/naturalearth/ne_110m_admin_0_countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 94 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

plot(st_geometry(countries), col="green")

Généralisation avec sf (ca marche mal)

countries2 <- st_simplify(x = countries, preserveTopology = TRUE, dTolerance = 2)

## Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
## st_simplify does not correctly simplify longitude/latitude data, dTolerance
## needs to be in decimal degrees

plot(st_geometry(countries2), col="green")

Généralisation avec rmapshaper (Youpie !!!)

countries2 <- ms_simplify(countries, keep = 0.1)
plot(st_geometry(countries2), col="green")

BIBLIO

BEGUIN Michèle, PUMAIN Denise, La représentation des données géographiques, Statistique et cartographie, coll. Cursus, Armand Colin, nouvelle édition 2000, 192p.
BERTIN Jacques, Sémiologie graphique, Monton-Gauthier-Villars, 1967, 1973, 432 p. Disponible à ce jour dans la collection Réimpression de l’EHESS, 1998.
BERTIN Jacques, La graphique et le traitement graphique de l’information, Flammarion, 1977, 280 p.
BRUNET Roger, La carte mode d’emploi, Fayard-Reclus, 1987, 270p.
ElementR (Groupe), R et espace, traitement de l’information géographique, 2014, Framabook.
LAMBERT Nicolas : Carnet de recherche neocarto
LAMBERT Nicolas, ZANIN Christine : « Manuel de cartographie. Principes, méthodes, applications », cursus, Armand Colin, 18 mai 2016, 224p.
Giraud Timothée : Carnet de recherche rgeomatic
ZANIN Christine, TREMELO Marie-Laure, Savoir faire une carte, Aide à la conception et à la réalisation d’une carte thématique univariée, coll. Sup géographie, Belin, 2003, 200p.

— R infos —

Version de R

R.version

##                _                           
## platform       x86_64-pc-linux-gnu         
## arch           x86_64                      
## os             linux-gnu                   
## system         x86_64, linux-gnu           
## status                                     
## major          3                           
## minor          4.4                         
## year           2018                        
## month          03                          
## day            15                          
## svn rev        74408                       
## language       R                           
## version.string R version 3.4.4 (2018-03-15)
## nickname       Someone to Lean On

Session

sessionInfo()

## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.1 LTS
## 
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] rmapshaper_0.4.0    osmdata_0.0.7.001   countrycode_1.00.0 
## [4] rnaturalearth_0.1.0 sp_1.3-1            cartogram_0.1.0    
## [7] cartography_2.1.2   readxl_1.1.0        sf_0.6-3           
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.18            lubridate_1.7.4        
##  [3] lattice_0.20-35         png_0.1-7              
##  [5] class_7.3-14            rprojroot_1.3-2        
##  [7] digest_0.6.17           V8_1.5                 
##  [9] mime_0.5                R6_2.2.2               
## [11] cellranger_1.1.0        plyr_1.8.4             
## [13] backports_1.1.2         rnaturalearthdata_0.1.0
## [15] evaluate_0.11           e1071_1.7-0            
## [17] geojson_0.2.0           httr_1.3.1             
## [19] highr_0.7               pillar_1.3.0           
## [21] rlang_0.2.2             lazyeval_0.2.1         
## [23] curl_3.2                geojsonlint_0.2.0      
## [25] rstudioapi_0.7          miniUI_0.1.1.1         
## [27] raster_2.6-7            geojsonio_0.6.0        
## [29] rmarkdown_1.10          rgdal_1.3-4            
## [31] foreign_0.8-69          jqr_1.0.0              
## [33] stringr_1.3.1           questionr_0.6.3        
## [35] shiny_1.1.0             rosm_0.2.2             
## [37] compiler_3.4.4          httpuv_1.4.5           
## [39] xfun_0.3                rgeos_0.3-28           
## [41] htmltools_0.3.6         tibble_1.4.2           
## [43] bookdown_0.7            jsonvalidate_1.0.0     
## [45] crayon_1.3.4            later_0.7.5            
## [47] grid_3.4.4              spData_0.2.9.4         
## [49] jsonlite_1.5            xtable_1.8-3           
## [51] DBI_1.0.0               magrittr_1.5           
## [53] units_0.6-1             rmdformats_0.3.3       
## [55] stringi_1.2.4           promises_1.0.1         
## [57] xml2_1.2.0              packcircles_0.3.3      
## [59] tools_3.4.4             abind_1.4-5            
## [61] yaml_2.2.0              maptools_0.9-4         
## [63] classInt_0.2-3          rvest_0.3.2            
## [65] knitr_1.20

GÉO-VISUALISATION AVEC RAction nationale de formationMATE SHS

GÉO-VISUALISATION AVEC RAction nationale de formationMATE SHS

Packages & données

Environnement de travail

Les packages

Packages dédiés à visualisation de données géographiques

Packages dédiés à la manipulation de données (spatiales ou non)

Packages fournisseurs de données

Téléchargez les données

Exo 1 - Bienvenue à Sète

Commandes de base

Introduction au SIG avec R / Manipuler les informations spatiales

A vous de jouer !

Exo 2 - Cartographier l’Occitanie

Objectifs

Chargement et préparation des données

Visualisation : carte en figuré proportionnel

Le package cartography

Éléments théoriques

Visualisation !

Visualisation : carte combinant figuré proportionnel et typologie

Visualisation : cartogramme de Dorling

A vous de jouer !

Exo 3 - Cartographier l’Europe

Objectifs

Chargement et préparation des données

Visualisation : cartes choroplèthes

A vous de jouer !

Visualisation : cartes carroyées

A vous de jouer !

Exo 4 - Cartographier le Monde

Objectifs :

Import du fond de carte

Import des données atributaires

Cartographie

Et si on ajoutait des lignes de discontinuités à cette carte?

BONUS

Utiliser des données OSM

Généraliser un fond de carte

BIBLIO

— R infos —

GÉO-VISUALISATION AVEC R

Action nationale de formation
MATE SHS

GÉO-VISUALISATION AVEC R

Action nationale de formation
MATE SHS